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Abstract 

We present analytic solutions to the spatially homogeneous ax- 
ion field equation, using a model potential which strongly resembles 
the standard anharmonic (1 — cos N9) potential, but contains only 
a piece-wise second order term. Our exactly soluble model for 9(t) 
spans the entire range [—ir/N,n/N]. In particular, we are able to 
confirm (i) Turner's numeric correction factors jjj to the adiabatic 
and harmonic analytic treatments of homogeneous axion oscillations, 
and (ii) Lyth's estimate || valid near the metastable misalignment 
angle ir/N at the peak of the potential. We compute the enhance- 
ment of axion density fluctuations that occurs when the axion mass 
becomes significant at T ~ 1 GeV. We find that the anharmonicity 
amplifies density fluctuations, but only significantly for relatively large 
initial misalignment angles. The enhancement factor is ~ (2,3,4,13) for 
in ~ (0.85, 0.90, 0.95, 0.99) x tt. 
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1 Introduction and model 



Axions are a consequence of the simplest solution to the strong CP problem 
P]. The complex Peccei-Quinn field ip has at temperatures above the QCD 
scale Aqcd an effective potential of the "Mexican hat" shape. Thus, tp = 
pe 1 /\/2 , with the classical value of p ~ the Peccei-Quinn breaking scale f a . 
The axion field a(x, t) corresponds to the angular variable / a 9 in the Peccei- 
Quinn field. Ignoring variations of p in the PQ field (which is reasonable 
for temperatures far below the PQ and inflation scales Q) we get the field 
equation 

□0(x,t) + V'(6(x,T(t))) = . 

The effective potential V for the ^-direction is flat for high temperatures, 
and the axion is massless at high temperatures. At low temperatures T ~ 
few x Aqcd ~ 1 GeV, the azimuthal part of the Mexican hat potential is 
tipped (N = 1) or rippled (N > 2) due to instanton effects ||, where A^ is 
the color anomaly of the PQ symmetry. The resulting curvature about the 
minimum is interpretable as the induced axion mass squared. The devia- 
tion of the potential from quadratic ^-dependence at larger 8 is the source 
of anharmonic effects. We will present a model that closely mimics the 
instanton-induced potential. Then we will analytically solve the model to 
obtain the time-evolution of fluctuations in the axion field 9, and ultimately, 
of fluctuations in the axion density. 

In principle, the initial axion perturbations could arise from stochastic 
values assigned to each acausal volume in the early, non-inflationary uni- 
verse, or from quantum mechanical field fluctuations within the causally 
connected horizon of the inflationary universe. In the acausal universe, one 
expects on average one domain wall per causal volume. These domain walls 
are disallowed by the observed smoothness of the universe's energy distri- 
bution, and so the inflationary origin of the axion density perturbations is 
preferred. However, the nature of the mechanism generating the initial fluc- 
tuations does not enter into our calculation to follow. This is because the 
appearance of the instanton-induced effective potential for the axion field 
occurs at T ~ 1 GeV, long after the epochs of inflation and PQ symmetry 
breaking. 

1 Models with changing p during inflation have been constructed ll. If they are true, 
they reopen the axion window to a wider cosmologically allowed range of axion masses. 
The calculations presented in this paper, however, are independent of the inflationary 
mechanism. 
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There is a subtlety to discuss. We will work with the homogeneous field 
equation, i. e. we will neglect gradient terms. A field inhomogeneity exists 
in a finite spatial region, and therefore the field has a nonzero gradient. The 
gradient will be negligible for small fluctuations, or for those fluctuations 
whose wavelengths are superhorizon during the time when the anharmonic 
effects come to dominate the expansion term, at T ~ 1 GeV. It is the 
evolution of these fluctuations that we consider. And it is these fluctuations 
that truly have no memory of their origin, e.g. whether they are entering 
(non-inflation) or re-entering (inflation) the horizon. Thus, as with all prior 
work, our analysis is strictly valid for small fluctuations satisfying 59 <C 
9. In fact, small scalar field fluctuations are expected in the inflationary 
scenario ||. Thus, axion fluctuations should be small but nonzero if the 
universe underwent an inflationary epoch that ended after T ~ / a and before 
T ~ Aqcd- Also, the potential domain wall problem which arises when 
9 = 8(1+59/9) approaches ±7r is mitigated in the small fluctuation universe. 
Recently, Kolb and Tkachev Q have initiated numerical integration of the 
partial differential equations which include the gradient terms; their results 
show that when gradients cannot be neglected, their effects may be dramatic, 
and that when fluctuations are large, domain walls will form. 

Perturbative calculations of the instanton-origin of the axion mass are 
only valid at temperatures T 3> Aqcd- In Turner's work a power law 



is fit to the perturbative instanton result for the potential's curvature at its 
minimum; values of n between 3.6 and 3.8 result. Since the perturbative 
treatment is not expected to be valid at just the temperature T ~ Aqcd 
where the mass is turning on, we will follow Turner and study a larger 
possible range for the power-law exponent n. Ultimately, the axion mass 
reaches its zero-temperature value and ceases to grow. 

9 parameterizes the different QCD vacua, whose energy densities, ne- 
glecting fermions^, are given by 



where So is the action of a QCD instanton, and the parameter K is given 



m a oc T 






by 



) 



2 



K = (m a f a /Nf = Urn. 



"7T 



l + z 



(79 MeV) 4 . 



(3) 
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The second equality in Eq. (||) is a current algebra result Q , relating the 
pseudo- Goldstone axion to the pseudo-Goldstone pion and to z, the ratio 
of the up quark to down quark mass. From it one finds, taking z ~ 0.5 
and N = 1, (m a /10 _5 eV) ~ (/ a /10 12 GeV) _1 . Astrophysical arguments 
[ 10 1 constrain the axion mass to lie in a rather narrow preferred "window", 



10 5 eV < m & < 10 3 eV. Further argument relates the present cosmic axion 
mass density (in units of the closure density f2 C rit) to the axion mass: £l a = 
/ \~ 7 / 6 

( uF^eV ) ' an ma = 10 eV axion closes the universe. The functional 
form of Eq. (0) is sensitive to the inclusion of fermions in the theory, although 



the periodicity of the potential is not 11 ] . Here we input a functional form 



that strongly resembles the familiar form of Eq. (||) , preserves the periodicity, 
but at the same time yields a homogeneous equation of motion which enables 
us to obtain an analytical solution for the field. We accomplish this by 
matching a parabola to an inverted parabola: 

J imJ(T)« 2 , |«| < f 

v(e) = \ to 2 ), |g|>j. < 4 > 

For convenience, we have assumed in Eq. (|j) that the color anomaly param- 
eter N in Eqs. (§) and (^) is one; our results extend trivially to the N > 2 
cases. Our model potential is compared to the cosine potential in Fig. 1. 
Previous analytic work on the evolution of the axion field was restricted 



to either the harmonic potential near the minimum at 9 = [12, 13|, or 



near the maximum at \9\ = -it [|2j]. In most cases, an adiabatic aproximation 

dm 
dt 



<C m 2 was imposed to facilitate a simple analytic solution O, 0] . Turner 



|l]] numerically integrated the equation of motion in the harmonic potential 
to test the adiabatic approximation, and in the anharmonic cosine potential 
of Eq. (|2|) for \9\ not too near n, to test the harmonic approximation. He 
found significant corrections. With our model potential, we are able to 
derive analytical solutions without succumbing to the harmonic or adiabatic 
approximations of prior work. Our solution thus improves and relates the 
prior work. 

We solve the field equations classically, which is valid during the epoch 



of perturbation growth, long after any possible inflation [14]. The only 
significant growth in density perturbations occurs at the time when the 
growing axion mass is comparable in magnitude to the diminishing Hubble 
parameter M. At earlier times the fluctuations are frozen by the overdamped 



equation of motion; and for later times, both classical [^, 13 1 and quantum 
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[15 1 calculations agree that fluctuation growth yields small overdensities. As 



we will show later in Fig. 4, the latter result derives from the fact that the 
late-time field has evolved into the harmonic region of the potential, where 
fluctuation growth does not occur. 

Once 9 is calculated, the energy density per physical volume is obtained 
from T^ v = fld^6d u 6 - g^C as p = T 00 = \fl9 2 + V(6), ignoring the 
gradient term. Eventually, the axion field evolves to small-amplitude, quasi- 
periodic oscillations in the harmonic region about the potential minimum. 
Then, to a good approximation, p is just V(0 O ut) ~ \ 7a m a$out> where 9 ou t 
is the oscillation amplitude and m a is the zero-temperature axion mass. 

If the initial fluctuations are sufficiently small, we can disregard spatial 
derivatives in the field equation and consider only the time-derivative part 
of the Klein-Gordon operator □. In a radiation-dominated FRW universe, 
this yields the usual equation of motion for the homogeneous "zero-mode" 
axion field (a(t) = f a 6(t)): 

m + Y t m+ml(T(t))9(t)=0, (5) 

if \9\ < 7r/2, and the same equation with 9 replaced by e = ±7r — 9 and m 2 
by — m 2 for \9\ > tt/2. Here the Hubble parameter, H(t), in the second term 
has been equated to l/2t, the value appropriate for a radiation-dominated 
FRW universe. 

Since the equation of motion is linear and homogeneous in 9 for \9\ < tt/2, 
the mean field 9 and the local field 9 = 9 + 59 evolve with the same time- 
dependence. Thus, the field contrast §=■ is time-independent, i.e. does not 

9 

evolve. The same conclusion follows for y at \9\ > tt/2. Thus, in our model, 
the evolutionary enhancement of y comes from the matching of solutions 

at 9 = ±7r/2. Growth occurs because 9 = 9 + 59 must be matched to ±7r/2 
at a time differing slightly from the matching time of 9. 

The argument of the mass in Eq. (||) requires a relation between time 
and temperature. In the radiation-dominated FRW universe, the two are 



inversely related by t oc T 2 . The relation is [16] t 1 = 3.3-y/g^j^, where 
Mp is the Planck mass ~ 1.2 x 10 19 GeV and g* is the effective number 
of degrees of freedom (with 7/8th's weighting for fermions) at temperature 
T. Counting photons, gluons, three neutrinos, two leptons, and three quark 
flavors, g* = 62 and t~ = 26-^. At T ~ 1 GeV, the value of 3H is 
~ 3 x 10~ 9 eV, which is much less than the zero-temperature axion mass. 
Thus, once the instanton-induced potential turns on at T ~ 1 GeV, its 
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curvature quickly dominates the Hubble term in the equation of motion. 

It is convenient to measure all time in units of t, the time (roughly) 
when the axion oscillations for the convex part of the potential first become 
undamped, t is defined implicitly by m a (i) = Jl, which yields T ~ 1 GeV 
and t ~ 10 -6 sec. In terms of this scaled, dimensionless time, the field 
equation becomes simply 

X(t) + j f x{t) ± ^t n X (t) = . (6) 

In the convex part of the effective potential, % denotes 9 and the plus sign 
holds for the last term; in the concave part of the potential, x denotes 
e = ±7r — 6 and the minus sign holds. The growth of the axion mass 
terminates when the mass attains its low temperature value m a (oo) at t c ~ 
(m a (oo)/m a (t)) 2 / n ; or equivalently, ~ (2/3m a 1 (oo)) 2 / n with t c and m~ 1 (oo) 
measured in units of t. 

If the initial value for \0\ exceeds vr/2, then each time the field amplitude 
passes through +tt/2 from above or below, there will be matching conditions 
to join together the and e solutions extracted from Eq. (||). The solutions 
to Eq. (|6|) are the Bessel and modified Bessel functions, respectively: 

( n + 2 s 
— (7) 
n + 2 j y ' 

, (8) 
n + 2 J ' y J 

2 The easiest way to see that scaling time to t is merely a convenience is to replace 
RW time in the equation of motion with conformal time r(t) = J -^fry (so that 
dr — dt/R(t)), and then to rescale the solution as \ = R( T )X- R( T ) is the scale factor of 
the conformal metric, which in the radiation dominated epoch equals R ^-- The resulting 
equation of motion then contains no explicit expansion term: 

'i±[R(r)m4T(r))] 2 x^0. 

Here, an overdot is a derivative with respect to conformal time. Since T scales as 1/R, we 
have, using Eq. (|), [R(t) m a (T '(r))] 2 = [7? m a (r )] 2 (^) 2n+2 , which makes it clear that 
any reference point chosen during the power-law growth period of the mass is equally 
valid. The solution to this e.o.m. is just x(r) = ^ Z ± _j_ ( R ° T °^( T °> (^)"+ 2 ), 
where Z is any Bessel function for the convex potential, and any modified Bessel func- 
tion for the concave potential. Connection with our RW-time solution is made by 
substituting the radiation dominance relation r/r = (t/to) 1 ^ 2 into \ — x/R to get 

X = {t /t) 1/i Z ± _i_ ^ fl ° T °'"° (To) ( 7 |-) 1 ^^. We see that our choice made in the text for 
the unit of time is just a particular value for the arbitrary reference point RoT 'm a ,(To) ■ 
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To match the initial condition at t <C 1, we will disregard the Bessel functions 
with negative order, since their singular behavior at t = implies that even 
a large component of this solution at the initial time t Q <C 1 quickly falls to 
insignificance by t = 1. However, whenever 6 evolves through ±7r/2 from 
above or below at a later stage, we have to match to the pair of independent 
solutions given above. The negative order solutions are the analogues of the 
singular mode arising (and customarily rejected) for the early epoch 

massless field discussed below. 

It is worthwhile at this point to examine the solution for the convex part 
of the potential, and to relate it to previous work done on homogeneous 
scalar field evolution. 

2 Mass turn— on with a harmonic potential 

The case of axion oscillations within the regime of a harmonic potential V oc 
6 2 has been discussed extensively in the literature Jl^, |l3[ . The following 
is a brief review of conclusions from these references. Before the axion 
oscillations become underdamped (i.e. before t), the mass term is negligible 
and the solutions to the field equation are simply 6 = constant, and 8 oc 
i -1 / 2 oc T. The constant solution is assumed to dominate when T cools 
to ~ lGeV. After the mass exceeds a few times the Hubble parameter, one 
takes the axion field to oscillate with frequency m(t), valid in the "adiabatic" 
approximation that m <C m 2 . 

With these approximations, the ansatz 

0(t) = A(t) cos(m(i)i) , (9) 

inserted into Eq. (||) yields A(t) oc i?~ 3//2 m -1 / 2 oc t~^ (R being the RW 
scale factor). With this adiabatic ansatz, the dependence of the axionic mass 
density on the initial misalignment angle is easily calculated. Apparently, 
mA 2 R? is an adiabatic invariant; call it J^// 2 . Since the axions are nonrela- 
tivistic, we see that the axion energy density is p a = ^m 2 f 2 A 2 = ^m&JA/R 3 
and the adiabatic invariant is just proportional to i? 3 /) a /m a , the axion num- 
ber per comoving volume. We also see that the energy density in the adi- 
abatic approximation scales as m a /i? 3 , and so the evolving energy density 
remembers the initial misalignment angle 9{ n = A(t < t) via the invariant 
J A = m a(T)9 2 n f 2 R(T)' i in the expression for the energy density: 

P,(T) = ^J A m a (T)/R 3 (T). (10) 
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The adiabatic approximation is invalid near t since rh a cannot be <C m 2 at 
the onset of oscillations. Furthermore, the adiabatic prediction of the ax- 
ionic amplitude evolution (time-independent J a oc R^p^/rria) implies that 
the axions behave like relativistic dust, with no axions produced or anni- 
hilated. This cannot be the case near or before t, not even for a purely 
harmonic potential, since the overdamped solution is <j) = constant in a 
physical volume, which means that axions are constantly produced in the 
comoving volume! On the other hand, the approximations are expected 
to rapidly become accurate at times shortly beyond t. To investigate the 
validity of the approximations made in the analytic work of [12[, Turner 



integrated Eq. @ numerically for a harmonic potential, as well as for 
the cosine-shaped anharmonic axion potential. For initial conditions within 
the harmonic part of the potential, he found a multiplicative factor f c {n) 
correcting the r.h.s. of Eq. fli~0|), which he fit with the linear expression 

f c (n) = 0.44 + 0.25n (numerical fit). 

(Recall that n is the power characterizing the axion mass turn-on in Eq. 
([!])). The correction factor is independent of #i n as a consequence of the 
homogeneous property of Eq. (||). With our solution to Eq. (||), we can in 
fact derive the correction factor. Let us do so. The solution for a harmonic 
potential with a power-law time-dependence (which is exactly our model 
potential Eq. (||), as long as |#i n | < tt/2), is 

1 / n + 2 \ 

/ 1 \ / 3 \"2^+4 1 (3t~\ 

0(t,|0 in | <7r/2)=0 in r 1 + — — ) — — t-±J_ 



i 



2n + 4 / \2n + 4j 2^+1 V n + 2 I 

The initial value 9{t = 0) = #i n determines the normalization factor. This 
result has, with different scaling of t, been published in ref. [13]. The pro- 



portionality constants are derived from the series expansion of the Bessel 
functions 



I v (z) j \2J ^n\T{n + u+l)\2 
The asymptotic form of Bessel functions of the first kind 

J±u(z) Z ^ J?- cos (z - (±2u + 1)^) + 0(z~ 3 / 2 ) (12) 

V 7TZ V 4/ 



3 When one inputs the perturbative instanton result n ~ 3.6 to 3.8, one obtains a 
correction factor f c (perturbative instanton) ~ 1.34 to 1.39, surprisingly close to unity. 
However, the validity of the perturbative instanton approximation, and therefore of the 
smallness of this correction factor, is suspect. 



8 



translates the J 1 solution into the late time expression 

2n+4 



2n + AJ \2n + AJ V 3vr 

( n + 2 \ 
3t~ TT 1 N \ 
n + 2 4 v n + 2 ; y v ; 

The n-dependent factors in the amplitude provide the improvement of the 
amplitude in Eq. (|)|), and the square of these factors correct the adiabatic 
axionic mass density presented on the r.h.s. of Eq. (JlQj) , viz. 

+ (U, 

Note that no approximations were made in deriving this non-adiabatic re- 
sult. There are two assumptions, namely, that the temperature dependence 
of the axion mass is given by Eq. ([[]) , and that the effective potential is purely 
quadratic (valid whenever the interactions, including self-interactions, of the 
scalar fields are sufficiently weak). 

The multiplicative correction factor f c is shown for < n < 10 in Fig. 2. 
The linear fit over the range of these values is given by f c = 0.456 + 0.247n, 
in very good agreement with the numerical results of [|l]] . The near-linearity 
of the correction factor in Eq. (|l4|) is evident in our figure over the entire 
range of n. 

That f c depends only on n and not on 0[ n implies that the probability 
distribution for 0-angles in a harmonic potential retains its shape for all 
times. That is to say, the relative over- and under-densities do not evolve. 
If the initial misalignment angle 9 = 8 + 59 is small enough to validate the 
harmonic approximation to the true axion effective potential, then we may 
write 

59 59 m 



for 



In particular, Gaussian amplitude fluctuations stay Gaussian after mass 
turn-on. This is not the case for fluctuations large enough to feel the an- 
harmonic effects. "Large enough" in our model potential given by Eq. (||) 
means \9\ > so that the field's initial value tastes the concave part of our 
potential. 
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3 Large initial misalignment angles: anharmonic 
evolution 



If |#i n | > we fit the initial condition 

e (t o <l) = e in (16) 

to the solution in @ with positive index, since the negative index solution is 
singular at t = 0. It may be that the initial value for e(t ) is also a stochastic 
parameter, smeared by the equal-time canonical commutator [9, 9] = i of 
the quantized theory. However, rather than complicate this analysis, we will 
assume that e(t Q ) is correctly given by omitting the negative order Bessel 
function at t Q . At the time t% when |eo(ii)| = § , eo must be matched to 
9i(t), which is a linear combination of the two functions in (Q); also, the 
first derivatives must be matched in the fashion eo(ti) = — For an 
initial condition not close to |#i n | = tt, matching at the field value ^ will 
only be necessary once, since the oscillation begins while the Hubble term 
is significant, acting as friction to reduce the field amplitude (as R~ 3 ^ 2 oc 
t -3 / 4 in the adiabatic approximation). The further evolution of the axion 
field then happens in the convex part of the potential, and is given by the 
functions in (^) . We will see that it suffices for quite a wide range of initial 
conditions to only calculate 9i(t). 

If, however, \9i n \ is close to tt, then the solution eo is metastable and so 
begins a late oscillation when the Hubble term is smaller; the matched solu- 
tion 9\ eventually evolves beyond the convex part of the effective potential. 
This necessitates a second matching, of 9\{t2) = — f to £2^2) and ^1(^2) to 
— ^2(^2)5 where €2(t) is a linear combination of functions in (||), and a third 
matching of £2(^3) = — § to #3(^3) and 62(^3) to —#3 (£3), where 9s(t) is a 
new linear combination of functions in ([?]). A schematic illustration of an 
oscillation pattern is given in Fig. 1. For large enough |0; n |, the harmonic 
minimum of the effective potential will be overshot arbitrarily often. This 
is because an initial field amplitude chosen near the maximum of the poten- 
tial will outlive the demise of the damping term in Eq. (0) which decreases 
rapidly with time. Said another way, a classical field amplitude of exactly tt 
will remain at tt and not oscillate at all! And so an initial amplitude slightly 
removed from tt will begin late oscillations when the damping/expansion 
is insignificant, and then oscillate for almost ever, in accord with energy 
conservation. 

For the construction of the individual solutions we use the notational 
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convention 

= vr - 9(t) > for t Q < t < h ; 

= 9(t) for ti < t < t 2 , with 6>i(ii) = e (ii 



eo(*) 



= c^e; ior ei <. e <. e 2 , wren t/i^i; = to^lj 

ea(*) = -7r-0(i) < forE 2 <t<t 3 , with e 2 (i 2 ) = #i(t 2 ) 

6 3 (t) = 9{t) for e 3 < t < t 4 , with 9 3 (t 3 ) = e 2 (t 3 ) 

= 7r — 6>(e) > for e 4 < i < t 5 , with e 4 (t 4 ) = ^3^4) 



e 4 (t) 



2 > 

7T 

2 

7T 

2 



(17) 

which is illustrated on the scale in Fig. 1. The series terminates at #at(e) = 
9{t) where 0jv(i) = (— l) - ~ f has no solution for t > eat. So the subscripts 
on e and 6 display a running count of the number of times t m that matching 
of solutions at the values 6 = ±ir/2 is required; equivalently, the index on 
9 m minus one is twice the number of times that the oscillating field has risen 
above the convex potential into the anharmonic regime. 
It is useful to introduce a new time scaling 



n + 2 

n + 2 



(18) 



In terms of z, the general forms of the solutions in eqns. (7) and (8) become 

9 m (z) = z- y [A m J u {z)+B m J^{z)] 
e n {z) = z- v [C n I u (z) + D n I_ v (z)) , 



(19) 



where v = (2n + 4) . Physical time t is obtained from z— time by inverting 
Eq. ©: 

'(n + 2)z^ (2Q) 



From Eq. (17) and the continuity of 9 it is clear that the connectivity con- 
ditions for the time derivatives are simply Q (with m an odd integer here), 



^m— 1 



(Z m ) 



i(z m ), and 6 m (z m+ i) 



(21) 



4 Since the common prefactor z~ v in Eq. ([[9]) is a C°° function at all z > 0, we may 
replace the e m and 9 m in the matching conditions with the simpler and equivalent z v e n — 
(C n Iv{z) + D„I- u (z)) and z v Q m — (A m J v (z) + B m J-v{z)). However, the recurrence 
relations on this page make it equally simple to work with the original functions. 
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To compute these derivatives, some recurrence relations are useful. For Z v 
any Bessel function, we have 

(z~"Z- v (z)) = z~ v Z_ {u+X) (z) , 

while 

±{z- v J v {z))=-z- v J v+1 {z) , 

and 

±(z-"I v {z))=z-"I v+1 (z). 

To compute the coefficients A m , B m , C n , D n , (m odd, n even), the pro- 
cedure is now clear. The initial condition eo(0) = e; n fixes the normalization 
of €q(z) to be 

eo(2f) = e in r(l + v) (J\ I v {z) ; 

n+2 

that is, Co = £i n r(l + v)2 v and Dq = 0. Next, Z\ = ^^t 1 2 is found by 

solving eo(z)\ z=Zi = ^. The conditions #1(211) = | and eo(-zi) = — #1(2:1) 
translate in terms of the notation in Eq. (|l9D into 



A x J v +\{z\)-B x J- v -\(z\) = CqI v+ x{z\) 

The inverse of 

M -- 



(22) 



is 



M" 



-i 



J V (Z) J-v{z) \ 

J v +x(z) -J-(v+q(z) J 

( J-(v+l)(z) J-u(z) 



W{J v ,J- v ;z) \ Jv+i{z) -J u {z) 



where W(J V , J- v ; z) = Jv(z)J_ {v+1) (z) + J^{z)J u+l {z) = is the 

Wronskian. Thus, 



M \ ixz\ I J_ {u+l) { Zl ) J- V (zi) \ ( \z\ 

Bi J 2sin7ri/ I J v +i{z\) -Ju(zi) ) \ C I v+ i(zi) 



(23) 



and 6*1(2;) = A±J u (z) +B\ J- V (z). Then we test whether or not the equation 
9\{z) = — f has solutions z > z\ . If it does not, then B\{z) is the final 
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solution for the field amplitude. If there are further solutions to the tran- 
scendental equation in z, then we name the smallest of these solutions Z2, 
and derive 



vrz 2 
2 sin 'kv 



(24) 



2 Z 2 



\ Iv+ifo) Iv{z 2 ) ) \ A l J u+ i{z 2 )- B 1 J_^ +1) {z 2 ) J ' 
using the Wronskian formula 

W(I u ,I- u ;z) = I v (z)Ii u+1) {z) - I_ v {z)I v +i(z) = . 

One can repeat this procedure as many times as needed to calculate all the 
further coefficients. For example, A3, B3, C4, and D4 are given by equations 
identical to ( |2~3| ) and p^ ) after some straightforward substitutions: 

(i) each subscript on A, B, C, D, and z is incremented by two, 

(ii) CqI u +\ in the column vector on the right-hand side of d23| ) is changed 
to C 2 I v +i + D 2 I-(y + \), and 

(iii) the matching value of the field, 5 in the column vectors on the right- 
hand sides of eqs. ( |23| ) and (^), is changed to the new matching value, — |. 

The final field amplitude Oj^(e- m ,t) = z~ u [A^J u + B^J„ U ] is generated 
in this recursive way. We show zi(ei n ) in Fig. 4. It is seen that the values 
of z\ hardly depend on n, whereas t\, which is an n-dependent reparame- 
terization of z, given in Eq. (^), depends sensitively on n. We also show 
the occurences and values of Z2 , Z3 , z^ and function of ej n , for values 

of n = 0, 3.7, 8. Several important inferences may be drawn from Fig. 4: 

(i) the convex region is overshot (for n = 3.7) when e in ~ 0.02 it, and double 
overshooting of the convex region occurs when e m ~ 10 ir; 

(ii) the z-time spent in the overshot concave region, Z3 — Z2, is relatively 
short; 

(iii) after a z-time equal to the largest Zi in the figure, the field oscillates in 
the harmonic region and no further growth of density perturbations takes 



5 Prior calculations used step-by-step integration of the field equation. When multi- 
oscillations are present, such an approach is computationally intensive. We proceed with 
a simple root-finder algorithm. We find qualitative agreement JIt]] with the earlier work 
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place; the largest Zi is relatively small on the cosmic time scale; 
(iv) we do not need to cut off the power-law growth of mass in our calcula- 
tions, since the anharmonic effects which govern density growth terminate 
before the axion mass reaches its low-temperature value: with the power- 

n + 2 

law growth of Eq. (|lj), the mass saturates at a z value of z c = (^4^)^ 2 = 
(7rp2")( m a(oo) /m a (t))~^~ . The exponent is always > 1. Thus, we ex- 
pect that z c > (^^)(m a (oo)/m a (t)). Since m a (oo) is > 10 -5 eV, and 

fn&ii) — 3H(T ~ lGeV) ~ 10~ 9 eV, we have z c ~ 10 4 , much larger than 
the occurrences of the matching times z\ through z$, even for infinitesimal 
values of e in . [] In practice, the basis {I U ,K U } offers an advantage over the 
basis {I± u } for expanding the e(z) solutions. This is because the asymptotic 
expressions for the former differ exponentially while those of the latter are 
in fact identical (being functions of v 2 rather than of v): 

and 

The use of K v rather than l— v gives much better stability for numerical 
work, and we have in fact used it. The substitution of K v for I_ u in the 
equations above for the A, B, C, D coefficients is effected with the following 
three changes: 

(i) I- v is replaced with K u ; 

(ii) l_(„_)_i) is replaced with —K u+ i, the relative minus sign arising from the 
fact that K y obeys the same derivative-recurrence relation as J u , opposite 
in sign to that of I— v \ 

(iii) the Wronskian W(I U , I—v'i z ) = — 2sm7TU is replaced with the Wronskian 
W(I u ,K u ;z) = I v {z)K u+1) {z) + K v {z)C +1 {z) = 1/z. 

The final values An and Bn are sufficient to determine the final density 
contrast, as we now show. From the asymptotic expansion of J v we can 

6 Once the mass has saturated at its zero-temperature value, the growth index n must 
go to zero, and the index and arguments of the Bessel solutions become fixed at 1/4 and 
|t, respectively. However, since the field evolution has entered its harmonic epoch prior 
to the onset of these solutions, these solutions are irrelevant to the growth of density 
perturbations. Moreover, since constant mass implies adiabaticity, there is no further 
increase in axion number. 
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compute the amplitude O ut at large z: 



6 N (z) = z~ v [A N J u (z) + B N J- V {z)\ 
= ont cos(z - 5), 



~ z- v J& [A N cos(z - f (1 + 2i/)) + B N cos(z - |(1 - 2z,))] 



with 



9 2 

7 out 



(27) 

(^tt) [(Ajv + B N ) 2 cos 2 (*£) + (A N - B N ) 2 sin 2 



2 / ra+2_\ n 
A 3 J 



n+3 

~+2 "+3 



2 



+ + 2AivB7V cos(^4) 



(28) 



and 



7T 

— + arctan 
4 



^4jv + Bn 2 



Note that the final coefficients An and Bn carry a memory of all the match- 
ing conditions, i.e. An = An{z\, Z2, Zn) and similarly for Bn] in turn, 
the Zi are predetermined by the value of the exponent n, and by the initial 
amplitude #i n = ±7r — t\ n . Finally, we may write down the asymptotic en- 
ergy density appropriate when the field is evolved into the harmonic region 
of the potential. It is just p a = |/ 2 "^ 2 (t)0 2 ut ) where # 2 ut is the squared field 
amplitude of Eq. (p8|), and m 2 (t) = m 2 (l)t n . 

It is interesting to compare our final solution in Eq. (|28| ) to the results 
of the harmonic approximation, with a time-independent axion mass, and 
with an adiabatically growing mass. The fixed mass solution to the har- 
monic Eq. (||) is simply 6 out = ^hJXfXif Jx (?n a f). The approximate 
answer is simple, but essential physics (mass turn-on, and anharmonicity) 
has been left out. The adiabatic solution to the harmonic potential, pre- 
sented in Eq. (|j) and explained thereafter, yields 9 2 oc (m a i? 3 ) -1 , and so 
# 2 ut = 9f n (i/t)~5~ . The asymptotic time-dependence of course agrees with 
our more careful answer, since in any model the field will evolve into the 
harmonic minimum of the potential at late times. However, the overall mag- 
nitude is different, because the harmonic approximation ignores important 
physics, at least for large 9 m . Finally, we comment that the solution to Eq. 
@ when the expansion term is omitted is y/iZ^. i (^y^ -5- )- with Z any 
Bessel (modified Bessel) function for the convex (concave) potential. The 
time dependence of this solution, and even the asymptotic behavior of the 
oscillation amplitude, ~ t _n//4 , show how much damping due to expansion 
affects our true solution. 
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The correction factor f c (0in) for the axion energy density from the adia- 
batic, harmonic solution is just the ratio of 9 2 ut in Eq. (28) to the harmonic 
result, 



7T 



n + 3 



s 2 f n + 2\ «+2 

Jc{"in 



A 2 N + B 2 N + 2A N B N cos( 



7T 



2n + 4' 



(29) 



This correction factor is plotted in Fig. 3, for three values of n. The correc- 
tion exceeds 2 when e^/n ~ 0.07, and exceeds 5 when €i n /ir ~ 3 x 10 _d . 

To quantify the final axion density contrast, = , we find it convenient 
to define its dependence on the initial field misalignment by 

Pa 2 0f n 9 in 

The density fluctuation 5p about the mean density ~p is defined by p = ~p+5p. 
The approximate expression in Eq. ( pCf ) is valid for small fluctuations. In the 
harmonic regime (9- m < ^ in our model), p a oc 9f n , and so £ = 2. For \9- m \ > 
t;, we expect anharmonic effects to retard the progression of the amplitude 
toward the minimum of the potential, thereby effecting a larger density 
contrast, £(6>i n ) > 2. In the extreme, we expect 5p a in our analysis to diverge 
as |#i n | — > 7T, because at the maximum of the potential there is no force to 
push or pull the field toward its minimum. However, in reality there are other 
limitations to this formalism at the potential maximum. For example, a field 
at a potential maximum will in general evolve in all directions characterized 
by a negative gradient. This then leads to domain walls. The energetics 
and evolution of domain walls are dominated by field gradient terms, and 
so are outside of our zero-mode formulation. As mentioned previously, the 
existence of walls conflicts with the deductions of observational cosmology. 

Since -Jr is just ■==, we have £(6> in ) = 2=^-70^. For small fluctuations, 



, „ ^v-m, - 2 

"out "out 



~ o ^ out _ o^ in ( d m ( flout ^ /o-n 

~ 9^ d9 in - 2 I-{ d ln(e in ) ) ■ {61) 

(Recall that = ±7r — 9- in .) Using our mostly analytic procedure we have 
computed values of £ for cases where the convex part of the effective potential 
is overshot once. We do so by computing 9 out as a function of e m within our 
model, and then taking the derivative numerically to determine £ as given 
in Eq. (||). 
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The advantages of our method over straightforward numerical integra- 
tion [jl], flTfl , are numerical stability and speed. Speed and stability are es- 
pecially advantageous for the most interesting case of a large initial mis- 
alignment angle. We use numerics only to find roots of a few transcendental 
equations and to perform 2 by 2 matrix multiplication to obtain the Bessel 



coefficients in Eqs. (23), (24), and their successors. After we obtain the val- 
ues of the final coefficients A jy and Bn, 9% ut is given by Eq. (|28|). We show 
this final enhancement of the density-fluctuations in Figs. 5 and 6. In Fig. 
5 the divergence arising from the factor e^T 1 in Eq. (|3l]) is scaled away and 
our results are compared with previous ones Q, and with our asymptotic 
solution, Eq. (|39|). In Fig. 6 we show the exact solution of our model for 
different values of n. The enhancement, £/2, is large for small e; n ; e.g. it is 
a factor of ~ (2,3,4,13) for e in ~ (0.15, 0.10, 0.05, 0.01) x vr. This is our most 
important result. 

If one were to abandon the small fluctuation condition 59 <C 9, then 
one might seek to compare the enhancements of densities arising from very 
different initial misalignment angles. The relative density enhancement is 
obtained by integrating Eq. (|30|). The result is 

— = exp / f(0) — 
Pi L-Mn,i 

From Fig. 6, it is evident that such an integral can become very large. Since 
grows rapidly with increasing large 9, we may write down an approxi- 

mate result: ^ « exp[£(6>2) ln(#2/#i)] = [ifA • However, large contrasts 
resulting from large initial field misalignments may not be trustworthy with- 
out inclusion of the gradient term JjJ . 

We will return to further discussion of the results displayed in Fig. 5 in 
the next sections. 



4 An asymptotic approximation 

We have convinced ourselves that for a large range of ei n values we get 
a reasonable approximation for the density contrast by including only the 
first occurence of 9(z), i.e. 0\{z), and by using the asymptotic expression 
for the Bessel solutions €q(z) and 9\{z) in the initial concave and the final 
convex regions, respectively. This may be so because for small ei n , z\ is 
3> 1, allowing the curvature of the potential to grow to a large value before 
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the onset of oscillations. This in turn implies a high oscillation frequency, 
so that the true solution simply does not spend much time in the concave 
potential (the duration of the €2 epoch is Z3 — Z2) after it overshoots the 
convex potential. These interpretations are supported by the zi, Z2, and Z3 
values in Fig. 4. Let us then pursue this approximate scheme of using only 
asymptotic eo and asymptotic 9\. 
Asymptotic eo is given by 

e (z) = C z- v I v {z) ^ C z-^=e z , (32) 



with Co = ei n 2^r(^ + 1). This equation is equivalent to the functional form 
obtained by Lyth 0], which he heuristically based on inverting the convex 
potential to a concave potential by replacing time with imaginary time, and 
applying the adiabatic approximation; equivalently, the asymptotic cosine 
solution is replaced by an asymptotic cosh. Q Thus, the satisfactory behavior 
of our asymptotic solution supports the approximations made in ||] . As long 
as the field remains in the concave regime, the field equation is linear and 
the density contrast remains constant: 

£ ^in 

From setting the asymptotic solution (^) equal to ^, we get the equation 
whose solution is the matching time z\\ 

°° zr ( " +l) e* = -• (33) 



2vr 1 2 



Even for this simple asymptotic solution (and even when the convex poten- 
tial is not overshot!), the late-time density fluctuations cannot be expressed 
as an analytic function of €i n , because the matching time zi(ei n ) is a solu- 
tion to the transcendental equation p^). Values for z\ obtained with these 
asymptotic formulae are indistinguishable from the solution to the exact 
equation ei n r(l + v){—) v 'I v (zi) = | shown in Fig. 4; from Eq. (25), the 
relative difference is (1 — 4^ 2 )/8zi. 

We may write the asymptotic solution in the convex regime as z~ v times 



a linear combination of the asymptotic forms of J± u given in Eq. (12), or 



7 The multiplicative correction factor to the adiabatically-evolving density fluctuation 
in the harmonic concave potential is exactly the same as the one obtained in Eq. (jlij) for 
the harmonic convex regime, which we found to be of order one. 
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alternatively, as 



6{z) ^1 A(e in , n)J -z~ {u+ ^ cos(z - 5) 



(34) 



with 5 = j + arctan( ^~ B * tan( i y : )). The asymptotic phase 5 is irrelevant 
for the calculation of the mass-density spectrum. Prom matching the field 
value in Eq. ( |34"| ) to 7r/2 at z = z\, we get 



A(e in ,n)J - cos(zi - (5) = ^Zi + 2 

V 7T Z 

To match the field derivatives at z\ we need 



(35) 



6{z] 



Z>1 



dz 

and from Eq. ( ^2| 



A(e in ,n)\j -z 

7T 



v + 



cos(z — 5) + sin(z — 5) 



(36) 



d i x Z»l 

^-e(^) — ► 
dz 



£ (z) . 



(37) 



From Eq. (pq) and the matching of the derivatives, Eqs. (^) and (p7|) with 
a relative minus sign, we eliminate the oscillating terms and obtain the time 
independent coefficient of the asymptotic amplitude 6: 



A(e in ,n) 



7T 



+ 2 (2i/+l) 2 
1- — 



Zl 



The n-dependence of A is implicit in v(n) and the ei n -dependence is implicit 
in z\{e- in ). The axionic energy density p a = \f 2 {0 2 + m 2 9 2 ) is then equal to 

, times an irrelevant time-dependent factor de- 



2 



2u+l 
1 



Au+2 , (2v+iy 



ducible from Eq. 



. Therefore, 



Pa(t) 



2v + l) 



1 + 



_2_ 



2(2^+1) 
— =2 



2 _ 4 1± 2 , (2i/+l) 2 

«1 Zl 2 



<5zi . , . 1 , c)zi 
= ~ (2j/ + 1)(1 + = ^ 

Zl Zl Zl 



for zi ^> 1. 

With the relation 



m 



5zi 
51 



^1 



+ * 



(38) 
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derived from Eq. ( p3[ ) , the relative density fluctuations in terms of the initial 
misalignment fluctuations are then given by 



Sp a (t) 



-(2i/ + 1) 

2/ 



f + 



1 

zT 



p a (t) z 1 -{v + 1 ) 
Notice that this means that for very small (implying large Z\) 

6p*(t) _ -{2v + l)5e- m 

pait) Zl ' 

or in the notation of eq. (|3C 



(39) 



Zl 7T 



We see then, that for 9 in near n, amplitude fluctuations evolve to give density 
fluctuations enhanced over the harmonic result by the factor 



2 7T — 9 m 2^1 



(40) 



In Fig. 5 this approximate result is presented and compared with the more 
exact treatment derived in the previous section, and with the previous esti- 
mate of Lyth Q . It is seen that this simple and powerful approximate result 
is about 30% low for ei n ~ 0.157T, and even better at smaller ei n . 



5 Discussion and conclusion 

To describe the whole history of the axion field and to compute the density 
fluctuations in the field, one begins with an assumed value for the mean 
misalignment angle #; n , and an assumed distribution about the mean for 
59i n , and then uses the formulation presented here to stretch the fluctuation 
spectrum during the time t ~ t. We have obtained a very simple expression 
for the evolution of axionic density fluctuations in the entire region £ 
[— 7T, 7r], valid for 59/6 <C 1 initially. By using a new model potential, we 
have been able to avoid the harmonic and adiabatic approximations, and 
still solve the field equations exactly. Our effective potential has the twin 
features of periodicity, and a power-law temperature dependence of the 
axion mass (with arbitrary power) as suggested by the finite temperature 
instanton approximation ||. 
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With our model we have extended the heuristic treatment of ref. J| away 
from 9 ~ ±7r. In ref. |0] an estimate of the oscillation period was used to 
estimate density fluctuations. Oscillation period, however, is an ill-defined 
concept around the time when m a ~ 3H and damping remains significant; 
and is ill-defined at any time for large initial angles where anharmonic effects 
are considerable. Nevertheless, good agreement between our result and the 
heuristic result is apparent in Fig. 5 (for the n = 4 power-law, which is an 
extrapolation of perturbative high-temperature instanton effects down to 
T). The anharmonicity amplifies density fluctuations, but only significantly 
so for relatively large initial misalignment angles. The enhancement factor 
we find is ~ (2,3,4,13) for 6 in ~ (0.85,0.90,0.95,0.99) x vr. 

There are possible complications to the axion evolutionary picture pre- 
sented here. These include 

(i) possible entropy production [p| and/or nonadiabatic axion mass evolu- 



tion |19|] during the QCD phase transition at T ~ Aqcd; 

(ii) possible damping of the axion oscillations due to interactions with other 

fields; 



(iii) possible axion "trapping" within quark or axion nuggets 21], in which 

the high-temperature phase survives. 

Let us briefly comment on each of these possibilities. 

The density fluctuations calculated in this paper are unaffected by en- 
tropy production or non-adiabatic mass change during the QCD phase tran- 
sition, provided the axion field is already in the harmonic regime by that 
time. This is because Eq. (|5|) stays linear for any temperature dependence 
of the axion mass, and the density spectrum is not distorted by evolution 
via a linear field equation, as discussed earlier. From the perturbative di- 
lute instanton approximation, Turner has extracted an estimate [U for the 
temperature at mass turn-on of 

/ f IN \ ~ a175 

r-o-s jp^r GeV - ( 41 ) 



V 10 12 GeV 



This temperature exceeds Aqcd ~ 100-300 MeV as long as / a /iV ~ 3 x 
10 14 GeV, or m a ~ 2 x 10~ 8 eV. The stronger bound / a < 10 12 GeV, or 
m a > 10~ 5 eV, that results from avoiding overclosure of the Universe with 
axions then provides a safety margin here, although models with hybrid 
inflation 0] or with considerable entropy generation after axion production 
22| allow for some relaxation of the overclosure bound. 

The interaction of the axion field with the thermal gluon field provides a 
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damping term [2C] in the equation of motion, Eq. (||). However, this damping 
is insignificant at mass turn-on and afterwards if f & /N ~ 10 10 GeV, or 
equivalently, m a < 0.6 x 10~ 3 eV. As mentioned in the introduction, axions 
more massive than 1CT 3 eV are disfavored by astrophysical constraints. If 
all uncertainties in the astrophysical lower bound on / a conspire, the / a — 
bound is lowered a bit to / a > 2 x 10 9 GeV @. In 

some models including 

non-standard astrophysics pffi , the bound relaxes to / a > 2 x 10 8 GeV. 
Thus it is possible that a 10 _3 -10 _1 eV axion might exist, in which case 
interaction-damping becomes important. 

In the scenario where axions become trapped in nuggets of high temper- 
ature phase, the low momentum axion modes which concern us here do not 
exist. The axion density fluctuation spectrum, basically flat in inflationary 
or stochastic models, is transformed into a k~ 2 spectrum within a limited k- 
interval determined by the horizon size at the QCD phase transition, much 
too small to affect structure formation [21]. Clearly, in this paper we must 
assume that axion trapping does not take place to any significant degree. 

The occurrence of large density fluctuations resulting from large initial 
field misalignments is possibly impacted by density-fluctuation bounds, e.g. 
COBE measurements. Evolution of the spectrum with the formalism pre- 
sented in this paper excludes regions in the m a -#i n plane. Since enhancement 
of density fluctuations is greatest at large 6, and since our results at large 9 
coincide with those of ref. ||] , our exclusions are the same as those of Lyth 
||. Given a model for the initial spectrum 0- m , or equivalently, for the initial 
field fluctuation spectrum =»-, these results may exclude particular models. 

Pin 

In summary, by solving the field equation exactly, we have provided a 
tool to evolve the amplitude fluctuations of the axion field. Our method 
is valid for any scalar field zero-mode (i.e. gradients are neglected) with 
temperature-dependent mass, provided background interactions may be ig- 
nored. The calculation was carried out for a radiation dominated RW metric. 
However, the generalization to curvature dominated, matter dominated, or 
other RW spacetimes is straightforward. Our evolution algorithm is ready 
and able to turn an initial fluctuation spectrum into a final density contrast. 
As the presently unknown initial fluctuation spectrum becomes more devel- 
oped theoretically, this algorithm will become more useful. Furthermore, 
recent work on the possible existence of axion miniclusters [25, 17, |?]] may 
find a complement in this work. 
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Figure Captions 

Figure 1: A comparison of our model potential (solid line) with the co- 
sine potential (dashed line); and a schematic drawing of the field oscillations 
starting from a large misalignment angle and evolving through the convex 
part of the potential, overshooting it once. The scale on the left illustrates 
the conventions we used to relate the functions e(z) and 9{z) to each other. 
The height of the potential is a function of time/temperature, but the shape 
is not. 

Figure 2: The multiplicative correction factor / c (n) given in Eq. ( |I~4T ) 
for a purely harmonic effective potential with mass turn-on according to Eq. 
(|l|). The dashed line is the linear fit which has been extracted numerically 
by Turner Q. 

Figure 3: The correction factor / c (#), for n = 0,3.7 and 8, to the 
axionic energy density of Eq. (|H]) predicted by the adiabatic and harmonic 
approximation . 

Figure 4: The "exact" z±, Z2, z%, 24, and Z5 from the transcendental 
equations with Bessel functions, vs. e^" for n = 0,3.7 and 8. The near 
equality of Z2 and 23 shows that the time spent in the overshoot solution 62 is 
small, especially since for z>l the z-axis represents a very stretched i-scale, 
suggesting that truncation of the solutions with 9\ is a good approximation. 

Note that the values of z\, Z2, etc. hardly depend on the choice of n, whereas 

„ 2 

tm = (^--2m)™+ 2 values obviously do. The higher is n, the earlier is the 
onset of oscillations, and the higher is the chance to overshoot the convex 
part of the effective axion potential. Using the asymptotic approximation 
to calculate Z\ gives a graph undistinguishable from this Figure. 
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Figure 5: £(Qm) n a Sin versus ei n , with n = 4. The enhancement factor 

for the density-contrast is The dot-dashed line is just the result in 
asymptotic approximation, given by Eq. (|39|), The solid line is the result 
using the exact solutions for the coefficients in Eqs. ( |23|) and fl24j) and in their 
successor equations. Here we considered the asymptotic solution only up to 
0\(z), but followed the exact numeric calculations up to 0%{z). Nevertheless, 
the agreement is fantastic, and shows that repeated reconnecting of the 
successive solutions is unnecessary. The discontinuity in the first derivative 
of the exact £(0- m ) a t £in ~ 0.06 is the result of the field overshooting the 
convex region and causing additional enhancement when e- m ~ 0.06. The 
dotted line is a previous estimate made by Lyth 0. 

Figure 6: £(0i n ) 

versus 6j n for various values of tl. We see that the 
enhancement depends much stronger on the initial misalignment angle than 
on the value of n. The values here need to be compared with £(0) = 2 in 
the harmonic approximation/regime. 
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